
% this function uses the output of parsePeriodicStimMovies_v4_1 where you
% get average gcamp responses to specific stimuli. You can average several
% sessions of this data by running parsePeriodicStimMovies_v4_1 and then
% registering the results to a common image (can be high resolution or low
% resolution). Processing will be with respect to THIS frame. This will the
% avGCAMP
%
% then need the positions of the vessels you're interested, again with
% respect to the widefield image everything is registered to. this is the 
% positions input. it should be two columns with first column the
% x-position and second column the y-position
%
% micronsPerPixel is again with respect to the image everything is
% registered to
%
% addition as of 12/19/2022:
% adding in input argument (fractionOfTotal which can be true or false) that will calculate instead of
% distance weighted activity, the ratio of distance weighted activity to
% total activity
%
% updates as of 2/22/2023
% replaced the distance calculation with a convolution between the gcamp
% image and the weighting kernel - should be much more computationally
% efficient 
function r = normalizedNeuralActivityDistance_v1_2(avGCAMP,positions,micronsPerPixel,scaling,activityClass,fractionOfTotal)

%% some presets for the function:
activityThresh = .2;
activityThresh2 = .5;           %this one is for looking for edge of activity
unweightedScaling = 20000;          %set a large number to calculate the "total" gcamp activity. Essentially creating a distance weight that weights everything in the window roughly equally
if nargin < 4 || isempty(scaling)
    scaling = 400;          %scaling factor in microns
end
if nargin < 5 || isempty(activityClass)
    activityClass = 'measured';
end
if nargin < 6 || isempty(fractionOfTotal)
    fractionOfTotal = false;
end
%% normalization step:
%need to fiddle with the cutoffs depending on the metric:
switch activityClass
    case 'simulated'
        %do nothing
    case 'measured'
        avGCAMP(avGCAMP<.20) = 0;
end

minval = zeros(size(avGCAMP,3),1);
for n = 1:size(avGCAMP,3)
    cur = avGCAMP(:,:,n);
    cur = cur(cur>0);
    if isempty(cur)
        minval(n) = 0;
    else
        minval(n) = min(cur);
    end
end
normGCAMP = (avGCAMP-median(minval));
normGCAMP(normGCAMP<0) = 0;
normGCAMP = normGCAMP/max(normGCAMP(:));
%% have to do morphological contraction to get rid of registration artifacts:
bw = normGCAMP>0;
SE = strel('diamond',3);
for n = 1:size(bw,3)
    bw(:,:,n) = imerode(bw(:,:,n),SE);
    bw(:,:,n) = imdilate(bw(:,:,n),SE);
end

normGCAMP(~bw) = 0;
%% calculate just the raw distance to the center of activity
dist2center = zeros(size(avGCAMP));
for n = 1:size(avGCAMP,3)
    [X,Y] = meshgrid(1:size(avGCAMP,2),1:size(avGCAMP,1));
    cur = avGCAMP(:,:,n);
    if max(cur(:)) == 0
        dist2center(:,:,n) = nan;
    else
        ind = find(cur == max(cur(:)));
        [y,x] = ind2sub(size(cur),ind(1));
        X = X-x; Y = Y-y;
        dist2center(:,:,n) = sqrt(X.^2 + Y.^2);
    end
end
r.dist2center = dist2center*micronsPerPixel;      %convert the distances in to microns

%% threshold the gcamp signal to exclude what is roughly noise:
thresholdedgcamp = normGCAMP;
thresholdedgcamp(isnan(thresholdedgcamp)) = 0;
thresholdedgcamp(thresholdedgcamp < activityThresh) = 0;
r.thresholdedgcamp = thresholdedgcamp;

%% get the distance-weighted gcamp activity for each position:
distanceWeightedActivity = nan(size(thresholdedgcamp,3),size(positions,1));

kernsz = ((((-2)*log(.001))^(1/4))*scaling)/micronsPerPixel;
[X,Y] = meshgrid((-kernsz:kernsz)*micronsPerPixel,(-kernsz:kernsz)*micronsPerPixel);
D = sqrt(X.^2 + Y.^2);
weightingkernel = exp(-.5*(D/scaling).^4);
indexedpositions = sub2ind([size(thresholdedgcamp,1),size(thresholdedgcamp,2)],round(positions(:,2)),round(positions(:,1)));
for n = 1:size(thresholdedgcamp,3)
    temp = conv2(thresholdedgcamp(:,:,n),weightingkernel,'same');
    distanceWeightedActivity(n,~isnan(indexedpositions)) = temp(indexedpositions(~isnan(indexedpositions)));
end

r.distanceWeightedActivity = distanceWeightedActivity;
 
%% get the distance to edge of clear activity:
thresholdedgcamp = normGCAMP;
thresholdedgcamp(isnan(thresholdedgcamp)) = 0;
thresholdedgcamp(thresholdedgcamp < activityThresh2) = 0;
r.thresholdedgcamp2 = thresholdedgcamp;
r.distance2activity = nan(size(thresholdedgcamp,3),size(positions,1));
indexedpositions = sub2ind([size(thresholdedgcamp,1),size(thresholdedgcamp,2)],round(positions(:,2)),round(positions(:,1)));
for n = 1:size(thresholdedgcamp,3)
    D = bwdist(thresholdedgcamp(:,:,n)>0);
    D = D*micronsPerPixel;
    r.distance2activity(n,~isnan(indexedpositions)) = D(indexedpositions(~isnan(indexedpositions)));
end

